rm(list = ls(all = TRUE))
library(foreign)
d <- read.dta("leaders.dta")
## Aggregate Estimates from MI

aggrMI <- function(mod){
    out <- list(NULL)
    beta <- matrix(NA, nrow = 5, ncol = 3)
    sigm <- matrix(NA, nrow = 5, ncol = 3)
    for(i in 1:5){
        beta[i, ] <- getME(mod[[i]], name = "beta")
        sigm[i, ] <- as.matrix(diag(vcov(mod[[i]])))
    }
    
    betabar <- round(colMeans(beta), digits = 3)
    outses <- rep(NA, 3)
    m <- 5
    for(i in 1:3){## See p. 53 in King et. al. 2001
        a <- sum((beta[, i] - betabar[i])^2) / (m - 1)
        b <- sum(sigm[, i]) / m
        c <- b + (a * (1 + (1 / m)))
        outses[i] <- sqrt(c)
    }
    out[[1]] <- betabar
    out[[2]] <- outses
    return(out)
}

## Table 1 Paper
library(texreg)
library(blme)
library(Amelia)
son1 <- glm(son ~ I(log(duration + 1)), family = binomial(link = "logit"), data = d[d$dep == 0, ])
son2 <- bglmer(son ~ I(log(duration + 1)) + (1 | country_id), family = binomial(link = "logit"), data = d[d$dep == 0, ])
son3 <- bglmer(son ~ I(log(duration + 1)) + (1 | country_id) + (1 | century), family = binomial(link = "logit"), data = d[d$dep == 0, ])
son4 <- bglmer(son ~ I(log(duration + 1)) + (1 | country_id) + (1 | century) + (1 | house_id), family = binomial(link = "logit"), data = d[d$dep == 0, ])
set.seed(444)
a.out <- amelia(d[, c("century", "country_id", "house_id", "name", "ascent_age", "ascent", "duration", "son", "brother", "next_deposed", "dep")], m =  5,
                bounds = matrix(c(5, 0, 87), nrow = 1, ncol = 3), ts = c("century"), cs = c("name"))
fA <- formula(son ~ I(log(duration + 1)) + I(log(ascent_age + 1)) + (1 | country_id) + (1 | century) + (1 | house_id))
son5 <- lapply(a.out$imputations, function(x) bglmer(fA, family = binomial(link="logit"), data = subset(x, x$dep == 0)))
son5aggr <- aggrMI(son5)

beta <- list(coef(son1), getME(son2, name = "beta"), getME(son3, name = "beta"), getME(son4, name = "beta"), son5aggr[[1]])
ses <- list(as.numeric(sqrt(diag(vcov(son1)))), sqrt(diag(vcov(son2))), sqrt(diag(vcov(son3))), sqrt(diag(vcov(son4))), son5aggr[[2]])
pval <- list(NULL)
for(i in 1:length(beta)){
    temp.beta <- beta[[i]]
    temp.ses <- ses[[i]]
    temp.pval <- 2 * (1 - pnorm(abs(temp.beta / temp.ses)))
    pval[[i]] <- as.numeric(temp.pval)
}

texreg(list(son1, son2, son3, son4, son5[[1]]), override.coef = beta, override.se = ses, override.pval = pval,
       stars = c(0.01, 0.05, 0.1), digits = 3, file = "table1.txt")

## Table 2 Paper
b1 <- glm(brother ~ I(log(duration + 1)), family = binomial(link = "logit"), data = subset(d, d$dep == 0))
b2 <- bglmer(brother ~ I(log(duration + 1)) + (1 | country_id), family = binomial(link = "logit"), data = subset(d, d$dep == 0))
b3 <- bglmer(brother ~ I(log(duration + 1)) + (1 | country_id) + (1 | century), family = binomial(link = "logit"), data = subset(d, d$dep == 0))
b4 <- bglmer(brother ~ I(log(duration + 1)) + (1 | country_id) + (1 | century) + (1 | house_id), family = binomial(link = "logit"), data = subset(d, d$dep == 0))
fA <- formula(brother ~ I(log(duration + 1)) + I(log(ascent_age + 1)) + (1 | country_id) + (1 | century) + (1 | house_id))
b5 <- lapply(a.out$imputations, function(x) bglmer(fA, family = binomial(link="logit"), data = subset(x, x$dep == 0)))
b5aggr <- aggrMI(b5)

beta <- list(coef(b1), getME(b2, name = "beta"), getME(b3, name = "beta"), getME(b4, name = "beta"), b5aggr[[1]])
ses <- list(as.numeric(sqrt(diag(vcov(b1)))), sqrt(diag(vcov(b2))), sqrt(diag(vcov(b3))), sqrt(diag(vcov(b4))), b5aggr[[2]])
pval <- list(NULL)
for(i in 1:length(beta)){
    temp.beta <- beta[[i]]
    temp.ses <- ses[[i]]
    temp.pval <- 2 * (1 - pnorm(abs(temp.beta / temp.ses)))
    pval[[i]] <- as.numeric(temp.pval)
}

texreg(list(b1, b2, b3, b4, b5[[1]]), override.coef = beta, override.se = ses, override.pval = pval,
       stars = c(0.01, 0.05, 0.1), digits = 3, file = "table2.txt")


## Table 3 Paper
nd1 <- glm(next_deposed ~ I(log(duration + 1)), family = binomial(link = "logit"), data = subset(d, d$dep == 0))
nd2 <- bglmer(next_deposed ~ I(log(duration + 1)) + (1 | country_id), family = binomial(link = "logit"), data = subset(d, d$dep == 0))
nd3 <- bglmer(next_deposed ~ I(log(duration + 1)) + (1 | country_id) + (1 | century), family = binomial(link = "logit"), data = subset(d, d$dep == 0))
nd4 <- bglmer(next_deposed ~ I(log(duration + 1)) + (1 | country_id) + (1 | century) + (1 | house_id), family = binomial(link = "logit"), data = subset(d, d$dep == 0))
fA <- formula(next_deposed ~ I(log(duration + 1)) + I(log(ascent_age + 1)) + (1 | country_id) + (1 | century) + (1 | house_id))
nd5 <- lapply(a.out$imputations, function(x) bglmer(fA, family = binomial(link="logit"), data = subset(x, x$dep == 0)))
nd5aggr <- aggrMI(nd5)

beta <- list(coef(nd1), getME(nd2, name = "beta"), getME(nd3, name = "beta"), getME(nd4, name = "beta"), nd5aggr[[1]])
ses <- list(as.numeric(sqrt(diag(vcov(nd1)))), sqrt(diag(vcov(nd2))), sqrt(diag(vcov(nd3))), sqrt(diag(vcov(nd4))), nd5aggr[[2]])
pval <- list(NULL)
for(i in 1:length(beta)){
    temp.beta <- beta[[i]]
    temp.ses <- ses[[i]]
    temp.pval <- 2 * (1 - pnorm(abs(temp.beta / temp.ses)))
    pval[[i]] <- as.numeric(temp.pval)
}

texreg(list(nd1, nd2, nd3, nd4, nd5[[1]]), override.coef = beta, override.se = ses, override.pval = pval,
       stars = c(0.01, 0.05, 0.1), digits = 3, file = "table3.txt")


save(son4, b4, nd4, file = "rdmodels_main.RData")
